We present a unified framework that integrates prior physical knowledge and kinematic constraints to advance data-driven simulation, system identification, and control in multibody dynamics applications. Departing from traditional methods that directly model system states using Fully-Connected Neural Networks (FCNN) or Recurrent Neural Networks (RNN), the proposed approach leverages the Neural Ordinary Differential Equation (NODE) paradigm [1] to model the derivatives of system states. This enables the handling of complex dynamics and index-3 differential-algebraic equations (DAEs) inherent in multibody dynamics. Referring to the figure, dropping for simplicity the subscript that stands for the time step, MBD-NODE models the independent accelerations via a neural net that is differentiable. Additionally, just like a physics-based simulator, NODEs can use the same set of parameters to evolve the system with states over time as a result of applying an input. In our approach, represents the set of independent generalized coordinates associated with the dynamic system, obtained via a coordinate partitioning step [2]. This ensures that physical laws and system constraints are inherently accounted for during simulation and control, without the need for separate inverse dynamics models or Lagrange multipliers. By integrating forward dynamics models with neural network derivatives, we can address inverse problems associated with system identification and model predictive control. The efficacy of our framework is demonstrated through numerical examples that showcase its ability to handle complex multibody dynamics systems with multiple. This unified approach provides an efficient and theoretically sound basis for simulation, system identification, and control applications in multibody dynamics. All models, data, and code supporting this work [3] are publicly available as open-source at https://github.com/uwsbel/sbel-reproducibility/tree/master/2024/MNODE-code.